#include <iostream>
#include <vector>
#include <cmath>
#include <fstream>
#include <string>
#include <cblas.h>
#include <lapacke.h>
#include "Spline.hpp" 
#include "UniqueSpline.hpp" 

// Example function to interpolate
double f(double x) {
    return 1.0 / (1.0 + x * x);  // f(x) = 1 / (1 + x^2)
}

// Subroutine to interpolate f using quadratic B-splines (degree 3)
void interpolateBSpline57(const std::vector<double>& x_vals, const std::vector<double>& y_vals) {
    BSpline spline(x_vals, y_vals, 3, "complete", -0.5, -20 / 101.0 / 101.0);
    std::ofstream file("spline_57.csv");
    if (!file) {
        std::cerr << "Error: Could not open file spline_57.csv for writing." << std::endl;
        return;
    }

    // Print out the interpolated values
    std::cout << "B-spline 57 Interpolation Results:" << std::endl;
    for (double x = x_vals[0]; x <= x_vals.back() + 1e-6; x += 0.1) { // Adjusted range
        double interpolated_value = spline.getSplineValue(x);
        std::cout << "x = " << x << ", Interpolated value = " << interpolated_value << std::endl;
        file << x << "," << interpolated_value << "\n";
    }
    file.close();
}

// Subroutine to interpolate f using quadratic B-splines (degree 2)
void interpolateBSpline58(const std::vector<double>& x_vals, const std::vector<double>& y_vals, double f1, double fN) {
    // Create a B-spline object (degree = 2)
    SpecialBSpline spline(x_vals, y_vals, f1, fN);
    std::ofstream file("spline_58.csv");
    if (!file) {
        std::cerr << "Error: Could not open file spline_58.csv for writing." << std::endl;
        return;
    }

    // Print out the interpolated values
    std::cout << "B-spline 58 Interpolation Results:" << std::endl;
    for (double x = x_vals[0]-0.5; x <= x_vals.back() + 1e-6+0.5; x += 0.1) { // Adjusted range
        double interpolated_value = spline.getSplineValue(x);
        std::cout << "x = " << x << ", Interpolated value = " << interpolated_value << std::endl;
        file << x << "," << interpolated_value << "\n";
    }
    file.close();
}

int main() {
    int N = 10;  // Number of data points
    std::vector<double> x1_vals(N);

    // Generate x1_vals (1, 2, ..., N)
    for (size_t i = 0; i < N; ++i) {
        x1_vals[i] = i + 1;
    }

    // Calculate the corresponding y values using the function f(x)
    std::vector<double> y1_vals(x1_vals.size());
    for (size_t i = 0; i < x1_vals.size(); ++i) {
        y1_vals[i] = f(x1_vals[i]);
    }

    // Interpolate using quadratic B-splines (degree 3)
    interpolateBSpline57(x1_vals, y1_vals);

    // Generate x2_vals (1.5, 2.5, ..., N-0.5)
    std::vector<double> x2_vals(N - 1);
    for (size_t i = 0; i < N - 1; ++i) {
        x2_vals[i] = i + 1 + 0.5;
    }

    double f1 = f(1);  // Boundary condition at the left
    double fN = f(N);  // Boundary condition at the right

    // Calculate the corresponding y2 values
    std::vector<double> y2_vals(x2_vals.size());
    for (size_t i = 0; i < x2_vals.size(); ++i) {
        y2_vals[i] = f(x2_vals[i]);
    }

    // Interpolate using quadratic B-splines (degree 2)
    interpolateBSpline58(x2_vals, y2_vals, f1, fN);

    return 0;
}

